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ABSTRACT: 

An algorithm for rational Chebyshev approximation based on computing 
the zeros of the error curve was investigated. At each iteration the pro- 
posed zeros are corrected by changing them toward the abscissa of the 
adjacent extreme of largest magnitude. The algorithm is formulated as a 
numerical solution of a certain system of ordinary differential equations. 
Convergence is obtained by showing the system is asymptotically stable at 
the zeros of the best approximation. With an adequate initial guess, the 
algorithm has never failed for functions which have a standard error curve. 
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1. Introduction 



We will consider the problem of computing the best approximation to a 



continuous function f (x) by a rational function R = a.x 1 / b .x^ 

mn h 1 j= 0 J 

where the measure of the error is the weighted sup norm, 



f - R I I = sup 
mn a<x<b 



f (x) - R (x) 
mn 

w(x) 



w(x) > 0 on [a,b] . A number of methods have been proposed, and a brief 
description and numerical comparison of some of them is given by Lee and 
Roberts [8] . 

We have implemented and analyzed an algorithm for computing the best 
rational approximation. This algorithm [7] has not yet appeared in the 
literature. The convergence rate is linear, and although the algorithm 
is relatively slow compared to some others, the algorithm was successful 
in instances where others, such as the Remes algorithm and Maehly’s second 
method, fail. For functions which have a standard error curve, the algorithm 
has never failed to converge to the best approximation, provided the initial 
approximant did not have a pole in [a,b]. 

The algorithm is formulated as the numerical solution of a system of 
m + n + 1 differential equations. The dependent variables are points of 
zero error in a proposed approximation. When there are m + n + 2 alternations 
in the error curve, the system of equations is at a rest point. Under the 
assumption of a relatively mild hypothesis, the algorithm is proved convergent 
by showing that the system of equations is asymptotically stable at that rest 
point . 
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2. A method for computation of R * 

mn 

We assume that the best approximation R * = P*/Q* exists, and is 

mn 

in reduced form. The basis for several algorithms, such as that of Remes 
(see, for example [3]) and Maehly’s second method [9] is the following 
characterization theorem, which can be found in may texts, e.g. Cheney [2] 



Theorem 1: In order that the irreducible rational function P/Q be a 

best approximation to f of the form R , it is necessary and sufficient 
that the error have a least 2 + max(m + deg Q, n + deg P) alternations. 

As with Maehly’s second algorithm, we assume that the error curve is 
"standard" in that it has exactly K = m + n + 1 points of zero error 
and K + 1 alternations. If the points of zero error are known, the 
best approximation can be computed as the rational function of the form 
R^ n which interpolates to the value of f at those points. 

Let z * < z * < ... < z * be the points of zero error for 
1 Z is. 



E* (x) = 



f(x) - R * (x) 
mn 

w(x) 



Let Z = (z. , . . . , z,,) be an approximation to Z* = (z z *) . Let 

IK IK 

R (x) be the rational function which interpolates to f at z , ...,z . 

mn l k 

Let z = a, z = b, and define N = sup | E (x) | , where E(x) = 

° K 1 Vr«\ 

R mn ^ . If Z is close to Z* , R exists and the N are 

y—r mn k 

w(x) 

finite. Let x^ be the point in suc h that |e(x^)| = , 

k = 1 K + 1 . For Z close to Z* these points must be unique. 
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Any of the variables with superscript asterisks will denote that variable 



for the best approximation. Note that N^^* - N^* = 0 , k = 1,2,...,K . 

Under our assumptions, if N^ + ^ - N^_ = ^ 9 ^ = 1>2,...,K > then (at least for 
Z in some neighborhood of Z*) Z = Z* , although in general this is not true. 

In a talk given at the joint SIAM-MAA meeting for the Northern California 
sections in February 1972, Dr. Milton W. Green of the Stanford Research 
Institute discussed an algorithm based on the above ideas [7] . Given an 
initial guess Z at Z* , one computes R , and then N , . . . ,N 



mn 



K+l # 



The value of each is then corrected by changing it so that its new 

value is nearer to or x ^. > as - N^ is positive or negative, 

respectively. That is, z^ is changed toward the point of largest (in 
magnitude) error in the interval z^^] * Green reported that he 

had had good success with the algorithm. 

In an attempt to systematize the method and to make it amenable to 
analysis for its convergence properties, we considered the basic idea in 
the following form. We formulated the method as a continuous (in the corrections 
to the z^) rather than a discrete problem. Consider the system of ordinary 
differential equations , 



= \ + i - \ 



k = 1, . . . ,K 



where z^ and N^ are as defined previously, and Z = Z^^ at t = 0 is 
an approximation to Z* , used as the initial condition. This system, when 
solved by Euler’s method, yields an algorithm similar to that proposed by 
Dr. Green. 
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The algorithm we study is based on a slightly different system of 
equations. Although the convergence properties are similar, we wished to 
remove the effect of linear transformations, and to incorporate some 
indication when the z^* bunch together, as happens when f has a large 
slope at some point. Consequently, we considered two somewhat modified 
systems of equations, neither of which seems to yield results markedly 
superior to the other. The first, and the one we analyze, is the system 

\+i ' N k 

(D z k = (x k+1 - x k ) , k = 1 K 

N 

where x is as defined previously, and N = max N , again with 
R l<k<K+l k 

Z = at t = 0 . The other system was 



( 2 ) 



where y. 



s 

( 



- N. 

. _ k+1 k | I i i ir 

z k _ I z k Y k l » k - 1,...,K 

N 



k i£ Vi- N k s0 • 



Vi 1£ N k+1 - N k * 0 > 



(0) 



and again Z = Z 
The factors x. 



at t = 0 . 

and 



k+1 X k ’ 



z^ - y^| are both an attempt to "slow" 



z^ when the z^* bunch together. In some cases (2) is superior to (1), 
and in others (1) is superior to (2) . We choose to analyze (1) because it 
seems to be more consistent in the optimum "time step" when solved by Euler* s 
method. The analysis of (2) is nearly the same. 

The point Z* is clearly a rest point of the system (1), and in the 
next section we give an analysis showing that under appropriate assumptions, 
(1) is asymptotically stable at Z = Z* . 
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There are many algorithms now possible, depending on the numerical method 
used to solve (1). Any method could be used, subject only to the appropriate 
choice of ’’time step”. However, we should bear in mind that the goal is not 
necessarily to solve (1) accurately, but rather to approach Z* closely. 

Thus, the use of Euler’s method for the solution. The choice of "time step” 
appears to be a matter of experience. It is desired to use a near optimal 
’’time step”, one which yields R * to the desired accuracy in a minimum 
number of time steps, or iterations. It is seen that the ’’time step” is a 
parameter similar to the parameter in the solution of elliptic boundary 
value problems by the alternating-direction implicit method [10]. It is 
doubtful that it can be taken arbitrarily large in our case, however. 

3 . Convergence 

The convergence of the algorithms possible in the setting of Section 2, 
when an appropriate ’’time step” is used, is determined by whether or not 
the system (1) is asymptotically stable at Z = Z* . We will make use of 
the following theorem, which is paraphrased slightly from the way it appears 
in most references, e.g. [1]. 

Theorem 2_: The system of equations 

(3) Z = A • (Z - Z*) + H(Z - Z*) , 

where A is a constant K x K matrix and H(Z - Z*) is a vector function 
which is small compared to Z - Z* , is uniformly asymptotically stable at the 
point Z = Z* if the eigenvalues of A all have negative real parts. 
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Thus, in order to analyze the system (1) we must put it in the form (3). 
Denote sign (E(x)) on (z^^, z^) by cr^ . By our assumptions about E(x) 
we can write 



K 

E (x) = G(x) n (x - z.) , 
i=l 1 



where G(x) is continuous and single signed if Z is close to Z* . Then 

K 

we have GCx^) [] (x^ -z.) , k=l,...,K+l . Recall also that 

i=l 

a = Z 0 < X 1 < Z l< *•• < X K < Z K < ^+1 = b * 



We now make an assumption about the dependence of G(x) on the z^ . 

As was done by Maehly and Witzgall [9], we assume that near the point Z* , 

the function G(x) does not depend very much on the z^ . (Note: In the 

£ 

case of approximating x by a polynomial of degree < K - 1 , with 
w(x) = 1 , we have G(x) = 1.) Then we have, for a given value of x , 



9E(x) _ x x E (x) 

dz. * - G(x) .n - *i> - - 

J 1=1 J 



Then, since x^ is the point of extreme error, even though x^ may change 



significantly with z^ , it is seen that the extreme value does not. Thus 



we have 



9Nj 

dz. 



\ E(x k> 
X k- Z j 



x. - z . 
k J 



since near Z* , N. ~ N. Then we 

k 



have 





4* H , k = 1 , . . • , K , 

IV 
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where H = 



is small compared to Z - Z* . Further simplification gives 



t\ 



\ H, 



K 



the expression 



z . - z 



(x * - x *V 
^ k+1 k ' 



j=l 



C x * - z *) (x * 
^ k+1 2 j ; kX k 



2.*) + H k ’ k 1 ’ 



. >K 



Thus, we have the system (1) in the form (3) with 

/ 

(4) A = 



<x k + i* - 2 j*> ( v - y> 



k,j !,•••, K 



Now we must investigate the location of the eigenvalues of the matrix 
A . To simplify the notation, we drop the asterisks from the variables. We 
will show that the matrix -A is a matrix of class |( (see Fiedler & Ptak 
[6]. These matrices are also known as M-matrices.) 

Let B be a square matrix with nonpositive elements, except possibly 
on the diagonal. Matrices of class K are a subset of such matrices, 
characterized by many equivalent properties [6]. The two equivalent pro- 
perties we shall use are given in 



Theorem 3 : The following properties of B are equivalent. (i) The real 
part of every eigenvalue of B is positive, ( i i ) Every leading principal 
minor of B is positive. 



Lemma 4 : The matrix A defined by (4) has negative elements on the diagonal 

and positive elements elsewhere. Further, the principal minors of -A , 
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are positive for f= 1,2,...,K , 



Det 



<x k+i - V 



<x k+l - 2 j> <x k ‘ 2 j)l 



> j !»•••»£ 



whenever x.. < z.. < x 9 < z„ < . . . < x < z < x 



Proof: We consider the element a 



K ^ K ^ K+l * 

(x k+i - V 



^ ( vi • y (x k - y 



The numerator 



is positive, and if j < k , then z_. < x_. + ^ < x^_ < x^ + ^ , and hence 

a^ >0. If j > k , we have x^ < x k+1 < x j < z j > and again a^_. > 0 . 

If i = k, x. = x. < z. <x.,_ = x, , hence a, . < 0 . We now consider 
k j 2 J+l k+1 

the principal minor of -A , 



Det 



(x kfi - y 



<x k+i - y <x k - y 



k, j - 



Removal of the factor (x - x ) from each row does not change the sign 

Kt X tC 



of the determinant, so we consider the 






Det 



( x k+i - y <x k - z j ) 



= D„ 



k, j = 1, . . . , £ 



The value of the latter determinant being nonzero is equivalent to the 



existence of a unique solution to the following interpolation problem: With 

1 



functions g^y).-— _ y) ^ . y) 



, k = !,...,£ , and given points 



(y . , w.) , j = 1,...,£ , find constants ot k = 1,...,£ such that 
3 2 x 

£ 



G „ W 



- E- 



O' 



'* ' ti (x k + i - V ( \ - V = * 
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Here we assume that none of the y. coincide with any of the x, , which 

3 k 

we have guaranteed in the case of interest* The above discussion says that 

is nonzero if the set of functions g^(y) , k = 1,...,£ is unisolvent on 

the permitted set of points. See Davis [5] for further discussion. 

The interpolation problem is known to be uniquely solvable, if and only 

if any linear combination of the g (y) , k = 1,...,£ which is zero at £ 

k £ 

distinct points is identically zero. Consider G^(y) = S a k^k^^ * Now 

k=l 



V y) ■ - 



£+1 

n 

k=l 



1 

(x k - y) 



£ 



E 

k=l 



£+1 

n 

j =1 

j*k+l 



(x j - y) = 



p t -i (y) 

Vi (y > 



> 



where ^ and are polynomials of degree £ - 1 and £ + 1 , res- 
pectively. We see that P , and hence can have at most £ - 1 dis- 
tinct zeros, unless is identically zero. Hence * 0 . 

Now let y , . . . ,y be variable, but such that each y. satisfies the 
1 £ 3 

condition stated for z. in the lemma, i.e., x. < y. < x . , - . Let y. , . . . ,y n 

J J J J+l 1 £ 

replace z^,...,z^ , respectively in , and note that x^ - y^ appears 

only on the diagonal. Hence by choosing y^ sufficiently close to x^ , the 

determinant is diagonally dominant, and is thus positive. By the continuity 



of as a function of z^ z^ , and the fact D f is never zero, we 

conclude that D > 0 . This completes the proof of the lemma. 



Application of Theorems 3 and 2 to the system (1) yields the fact that 
(1) is asymptotically stable at Z = Z* . 
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4. Numerical Implementation, 

A version of the algorithm, which we will call "Algorithm G", was 
implemented in double precision Fortran on the IBM 360/67 at the Naval 
Postgraduate School. We used Euler’s Method to solve (1) numerically, using 
the initial guess to be the zeros of the Chebyshev polynomial of 

degree K , translated to the interval [a,b]. The rational function R 

mn 



was found by solving the linear system obtained by requiring R to 

mn 

interpolate to f at the . The IMSL^ routine LEQTIF was used to 

solve the system. 

The extreme values N^ were found by the method suggested by Maehly 

and Witzgall [9]. A search is made for a "turning point" using the 

previous value of as an initial estimate, using steps of 

h = . 015(z. - z. -). When three points have been computed so that 
k k k-I 

|E(x) | is largest at the middle point, the value of x^ is approximated 
by passing a parabola through the three points, and finding its extreme 
point. The alternative is to convert the problem to a discrete one by 
evaluating the error at a fixed number of points, as did Lee and Roberts [8]. 
We feel our method is probably faster, especially in the latter stages, 
although it assumes the error curve has but one local extreme in each 
interval (z^ z^) . We believe the method to be more accurate for smooth 

problems, as well as preserving the continuity of the original problem. 



1 IMSL - International Mathematical and Statistical Libraries, Inc., 
6200 Hillcroft, Suite 510, Houston, Texas 77036 
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Algorithm G was tested by running a variety of problems. The same 

2 

problems were also run using the Remes algorithm and Maehly's second 
method; the latter program was written by the author. The problems ranged 
from M easy M , such as exp(x), T(x) , and log x, to "hard" problems such 



function w(x) = v/x - for relative fits. In addition, the function r(x) = 



This latter function is due to Rutishauser (see Cody & Stoer [4], p. 179) 
and is difficult because the initial guesses usually used (data from 
appropriate Chebyshev polynomial) with the Remes algorithm lead to 
approximants with a pole in [-1,1] . 

The approximations of the form R^, ^13’ ^22* anc * ^42 were usually 

attempted, although in some specific instances others were computed. In 
the case of the "easy” functions all three algorithms worked well, with 
the Remes algorithm converging very rapidly, of course. The Remes algorithm 
and Maehley's method failed on one or more of the M hard" problems, while 
Algorithm G, for appropriate "time step", did not fail. We note that one 
iteration of the Remes algorithm requires the solution of a non-linear 

system, Maehly’s method requires the solution of two linear systems, while 
Algorithm G requires the solution of one linear system. 

The difficulties with Rutishauser f s function r(x) were investigated 
more closely. Aside from the problem of poles, the function has a large 
slope near x = 0 , which apparently causes difficulties. The programs 



2 The IMSL routine IRATCU was used. This is a Fortran version of procedure 
Ck&byAhdV due to Cody, Fraser, and Hart [3]. 



as v/T on [0,1] and \fx on [-^, 1 ] and [j^-,1] , the latter two with weight 



(arctan 8x) 




with w(x) = 1 , on [-1,1] was attempted. 
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for the Remes algorithm and Maehly's method were modified to accept input 
initial guesses. For initial guesses at the extremes which were accurate 
to seven significant digits, the Remes algorithm failed. For initial guesses 
at the interpolation points which were accurate to seven significant digits, 
Maehly's method failed. No particular difficulties were experienced by 
Algorithm G. 

With regard to the possibility of poles in the initial approximant, 
we have discovered that while theoretically the method fails, the numerical 
algorithm may be able to recover. For example, when approximating sin x 
on [0, 4.1] with = Pq/Q^ > using the Chebyshev points as initial 

guesses gives an initial approximant with a pole near x = 1.7. However, 
because of the approximation to the extremes, the routine recovers, forces 
the pole out of the approximation, and then converges to the correct result. 
This may not be a general rule, however, but is an interesting example of 
how robust the algorithm can be. 

Having satisfied ourselves that the algorithm works quite well, a study 
was made of how one should choose the "time step". We found that "time step" 
At = .20 was usually (not always) small enough for convergence. The optimum 
value for At is dependent on the problem, and sometimes is significantly 
larger than .20. As might be expected, the optimum value is larger than 
required for accurate solution of (1), and it is better to underestimate 
At than to overestimate it. Table 1 gives the number of time steps (or 
iterations) required for convergence of various approximations to exp(x) 
on [0,1], on [0,1], and r(x) on [-1,1], versus At . The iterations 

were stopped when successive approximations were obtained whose respective 
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coefficients had relative differences of less than 10 

Because the solution of (1) is close to the solution of Z = A(Z - Z*) 
for large times, the convergence rate is seen to be linear. Further, we 
can see in Table 1 that convergence is slow even for smooth functions 
such as exp(x), where 10-20 iterations are required. In the case of \/x , 
where most of the z^ tend to bunch near zero, significantly more iterations 
are required, the number increasing with K . 



Table 1: Search for optimum At 



exp (x) \/x~ r (x) 



At 


V^I 


V^2 


P /Q 
r v i 


P 2«2 


V Q 1 


Vq 2 


v q 2 


.10 


36 




64 










.125 






51 


>100 


29 


64 


61 


.15 


23 


47 


41 


97 


22 


49 


53 


.175 






44 


82 


17 


39 


42 


.20 


16 


35 


>100 


>100 


22 


failed 


failed 


.225 


14 








32 







.25 12 28 

.275 19 

.30 36 23 

.325 

.35 19 

.375 18 

.40 16 

.425 16 
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5. Conclusion. 



Algorithm G appears to be very robust. When coupled with the appropriate 
"time step", it is likely infallible, provided the initial guess did not 
yield an approximant with a pole in the interval. Even then, the version 
of the algorithm we implemented has recovered in specific instances. 

While this algorithm cannot compete with the Remes algorithm on the 
basis of speed, our tests show its speed to be comparable to Maehly’s 
algorithm, which shows up quite well, on that basis, in the Lee and Roberts 
study. One suspects the Lee and Roberts timing is biased since cases where 
the algorithm was successful were likely to be the relatively easier (and 
faster) cases. 

The simplicity of the algorithm coupled with its high success rate, 
make it worthy of consideration for computing approximations which result 
in failure of the Remes algorithm. The principal disadvantage seems to 
be the assumption that there are no more than m + n + 1 zeros of the 
error curve. Degenerate cases can be handled by increasing the degree of 
the numerator and/or denominator so that the error curve becomes standard. 

If the cause for more than m + n + 1 zeros is not degeneracy, the algorithm 
may fail. 

6 . Acknowledgement 

The author wishes to acknowledge helpful discussions with Professors 
A. L. Schoenstadt and D. A. Archer. 



14 



REFERENCES 



1. S. Barnett and C. Storey, Matrix methods in stability theorey , 

Barnes and Noble, Inc., New York, 1970. 

2. E. W. Cheney, Introduction to approximation theory , McGraw-Hill 
Book Co., New York, 1966 

3. W. J. Cody, W. Fraser and J. F. Hart, n Rational Chebyshev approximation 
using linear equations 11 , Numerische Mathematik 12(1968) 242-251. 

4. W. J. Cody and J. Stoer, n Rational Chebyshev approximation using inter- 
polation", Numerische Mathematik 9(1966) 177-188. 

5. Philip J. Davis, Interpolation and approximation , Blaisdell Publishing 
Co., Waltham, Mass., 1963. 

6. Miroslav Fiedler and Vlastimil Ptak, "On matrices with non-positive 
off-diagonal elements and positive principal minors" Czechoslovak 
Mathematical Journal 12(1962) 382-400 (English). 

7. Milton W. Green, "A new algorithm for rational approximation" delivered 
at joint meeting of the Northern California sections of SIAM and MAA, 
February 5, 1972. 

8. C. M. Lee and F. D. K. Roberts, "A comparison of algorithms for rational 

approximation", Math. Comput. 27(1973) 111-121. 

9. Hans J. Maehly and Christoph Witzgall, "Methods for fitting rational 
approximations, Parts II and III", Journal ACM 10(1963) 257-277. 

10. D. W. Peaceman and H. H. Rachford, "The numerical solution of parabolic 
and elliptic differential equations", J. SIAM 3(1955) 28-41. 



15 



Distribution List 



No. of copies 



Defense Documentation Center 12 

Cameron Station 
Alexandria, Virginia 22314 

Library 2 

Naval Postgraduate School 
Monterey, California 93940 

Dean of Research Administration 2 

Naval Postgraduate School 
Monterey, California 93940 

Professor Richard Franke 5 



Department of Mathematics 
Naval Postgraduate School 
Monterey, California 93940 

Professor Craig Comstock 5 

Department of Mathematics 
Naval Postgraduate School 
Monterey, California 93940 

Professor David A. Ault 1 

Department of Computer Science 

Virginia Polytechnic Institute 

and State University 

Blacksburg, Virginia 24061 

Professor W. Max Woods 1 

Chairman, Department of Mathematics 
Naval Postgraduate School 
Monterey, California 93940 

Professor R. E. Barnhill 1 

Department of Mathematics 

University of Utah 

Salt Lake City, Utah 84112 

Professor A. L. Schoenstadt 1 

Department of Mathematics 
Naval Postgraduate School 
Monterey, California 93940 

Professor D. A. Archer 1 

Department of Mathematics 
Naval Postgraduate School 
Monterey, California 93940 



17 



No. 



Dr. Milton W. Green 
Stanford Research Institute 
Stanford, California 94305 

Dr. Richard Lau 

Office of Naval Research 

Pasadena, California 91100 

Dr. Leila Bram 

Director, Mathematics Program 
Office of Naval Research 
Arlington, Virginia 22217 

Professor Frank Lether 
Department of Mathematics 
University of Georgia 
Athens, Georgia 30601 

Professor 0. G. Johnson 
Department of Computer Science 
University of Houston 
Houston, Texas 77004 

Mr. Duane Stevenson 
Kaman Nuclear 

1700 Garden of the Gods Road 
Colorado Springs, Colorado 80907 

Mr. W. J. Dejka, Code 4300 

Naval Electronics Laboratory Center 

San Diego, California 92152 

Professor S. E. Weinstein 
Department of Mathematics 
Arizona State University 
Tempe, Arizona 85281 

Professor D. E. Olson 

Department of Electrical Engineering 

University of Idaho 

Moscow, Idaho 83843 

International Mathematical 
and Statistical Libraries, Inc. 

6200 Hillcroft, Suite 510 
Houston, Texas 77036 



of copies 
1 

1 

1 

1 

1 

1 

2 

1 

1 



18 



Unclassified 



SECURITY CLASSIFICATION OF THIS PAGE (Whan Date Entered) 



REPORT DOCUMENTATION PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1 . report number 

NPS-53FE73111A 


2. GOVT ACCESSION NO. 


3. RECIPIENT'S CATALOG NUMBER 


4 TITLE (and Subtitle) 

On the Convergence of an Algorithm for Rational 
Chebyshev Approximation 


5. TYPE OF REPORT & PERIOD COVERED 

Technical Report, 1973 


6. PERFORMING ORG. REPORT NUMBER 


7 . author^; 

Richard Franke 


8. CONTRACT OR GRANT NUMBERS 


9 PERFORMING ORGANIZATION NAME AND ADDRESS 

Foundation Research Program 

Naval Postgraduate School, Monterey, CA 9394C 


10. PROGRAM ELEMENT. PROJECT. TASK 
AREA A WORK UNIT NUMBERS 


M. CONTROLLING OFFICE NAME AND ADDRESS 


12 report DATE 

November 1973 


13. NUMBER of pages 

20 


U. MONITORING AGENCY NAME A ADDRESS^/ different from Controiilng Office) 


15 SECURITY CLASS, (of thla report) 

Unclassified 


15a. decl assi fi cation/ downgrading 
schedule 


16. DISTRIBUTION STATEMENT (ot tht» Report) 

Approved for public release; distribution unlimited. 


17. DISTRIBUTION STATEMENT (of the ebetract entered in Block 20, If different from Report) 


10 . supplementary NOTES 


19. KEY WORDS (Continue on reveree aide If necaaaary mnd identify by block number) 

Rational Approximation 
Chebyshev Approximation 
Stability 
M-Matrices 


20 . ABSTRACT (Continue on raverae aide It naceaamry mnd identity by biock number) 

An algorithm for rational Chebyshev approximation based on computing 
the zeros of the error curve was investigated. At each iteration the 
proposed zeros are corrected by changing them toward the abscissa of the 
adjacent extreme of largest magnitude. The algorithm is formulated as a 
numerical solution of a certain system of ordinary differential equations. 
Convergence is obtained by showing the system is asymptotically stable at 
the zeros of the best approximation. With an adequate initial guess, the 
algorithm has never failed for functions which have a standard, error curve. 
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